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We present a reduced model for the transition to turbulence in shear flow that is 
simple enough to admit a thorough numerical investigation while allowing spatio- 
temporal dynamics that are substantially more complex than those allowed in pre- 
^^ I vious modal truncations. 

Cn . Our model allows a comparison of the dynamics resulting from initial perturba- 

tions that are localised in the spanwise direction with those resulting from sinusoidal 
perturbations. For spanwise-localised initial conditions the subcritical transition to 
a 'turbulent' state (i) takes place more abruptly, with a boundary between laminar 
and 'turbulent' flow that is appears to be much less 'structured' and (ii) results in 
a spatiotemporally chaotic regime within which the lifetimes of spatiotemporally 
complicated transients are longer, and are even more sensitive to initial conditions. 
The minimum initial energy Eq required for a spanwise-localised initial pertur- 
bation to excite a chaotic transient has a power-law scaling with Reynolds number 
^ Eq ^ Re^ with p w —4.3. The exponent p depends only weakly on the width of 

the localised perturbation and is lower than that commonly observed in previous 
t^ ' low-dimensional models where typically p « —2. 

. ^ . The distributions of lifetimes of chaotic transients at fixed Reynolds number are 

£^ ' found to be consistent with exponential distributions. 

r^ ' Keywords: fluid flow; turbulence; dynamical systems 

1. Introduction 

^ I The transition from laminar to turbulent states has been a central problem in fluid 

^~r: ' mechanics for many decades. Since the 1960s, promising lines of attack have been 

jy~ I opened up through the use of ideas from dynamical systems theory coupled to the 

f^ ' thorough investigation of reduced versions of the Navier-Stokes equations. Per- 

I ' I haps the best-known example of such a reduction is the derivation of the Lorenz 

f— ^ ■ equations as a model of the onset of thermal convection in a layer of fluid heated 

from below (Lorenz 1963). Although the set of three nonlinear ODEs that comprise 
the Lorenz 1963 model displays a wealth of interesting dynamical behaviour, little 
of this is relevant to the original fluid-mechanical problem. However, similar ap- 
proaches yield convincing agreement over wide ranges of parameter values in other 
fluid-mechanical situations, for example thermal convection in the presence of a 
j^ i magnetic field, and the onset of Taylor vortices in the flow between rotating coax- 

ial concentric cylinders. In situations such as these the flow undergoes a series of 
bifurcations before becoming chaotic, in a sense that can be given a clear meaning 
in terms of the behaviour of the reduced model comprising nonlinear ODEs. 
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In contrast, the transition to turbulence in shear flows appears not to proceed 
through a sequence of bifurcations, but to be linked to the appearance, at a critical 
Reynolds number Rtc, of a chaotic saddle in phase space: a comphcated collection of 
unstable equilibria and time-periodic orbits which results in ever longer transients 
before the flow relaxes to a purely laminar state. For a review, and substantial 
numbers of references to the literature (although with an emphasis on pipe flow), 
see Kerswell (2005). The significance of Rcc is that for Re < Rcc the laminar 
state is a global attractor and trajectories appear to evolve rapidly towards it, 
while for Re > Rbc other (usually unstable) invariant sets exist in phase space. 
These new invariant sets cause increasingly long transient excursions to take place 
before relaminarisation occurs, for initial conditions that are sufficiently far from 
the laminar profile. Small perturbations to the laminar state still decay rapidly 
towards it, and there appears to be a distinct boundary separating the behaviour 
of trajectories into those that relax rapidly and those which undergo long transient 
excursions. This laminar-turbulent boundary is sometimes referred to as the 'edge of 
chaos'. Despite these differences in phenomenology and dynamics, reduced models 
constructed along the same lines as the Lorenz model have provided substantial 
insight into both the physical origin of this transition to self-sustaining complicated 
motion, and the mathematical organisation of equilibria and periodic orbits inside 
the chaotic saddle. 

The study by Waleffe (1997) provides the direct inspiration for the present work. 
Waleffe showed that a low-order reduced model could be constructed that elucidated 
the different elements of a self-sustaining process (SSP) that allowed sufficiently 
large deviations from the laminar flow profile to persist indefinitely. The SSP can be 
briefly described as follows. Weak streamwise vortices (i.e. vortical rolls whose axes 
are aligned with the primary flow direction) distort the streamwise velocity profile 
by moving high and low-speed fluid around. This distortion generates streamwise 
streaks of fluid that are moving faster and slower than the fluid around them. The 
streaks are unstable to modes which create vortical eddies oriented in the wall- 
normal direction, orthogonal to the streamwise vortices, and the resulting three- 
dimensional re-organisation of the flow reinforces the streamwise vortices. 

Waleffc's model simplified the Navier-Stokes equations first by considering, not 
plane Couette flow between rigid boundaries, but a modified problem in which the 
boundaries are stress-free and the laminar profile is sinusoidal, sustained by an 
artificially-applied pressure term. It appears that the physics of the SSP is rather 
insensitive to these modifications. The second set of simplifications concern the 
Galerkin expansion of the velocity field in all three directions: streamwise (x), wall 
normal (y) and spanwise (z). Thus Waleffe considered a model comprising 8 of the 
lowest-wavenumber Fourier modes in this Galerkin expansion. These 8 modes were 
chosen in order to capture the central elements of the SSP, and to be self-consistent 
in the sense that the nonlinear (quadratic) interactions between these 8 modes 
preserve energy just as the full advective nonlinearity in the Navier-Stokes equa- 
tions does. Having projected out the spatial dependence of the dynamics onto these 
modes, the problem reduces to a far simpler set of ordinary differential equations 
(ODEs) for the time-dependent mode amplitudes. Subsequent work, in particular 
by Eckhardt & Mersmann (1999) and Mochlis, Faisst and Eckhardt (2004, 2005) 
extended Waleffe's model to include an additional physical insight: that the basic 
laminar profile of the shear flow will itself be modified by the nonlinear interac- 



Article submitted to Royal Society 



Transition in a ID shear flow model 



tions with streamwise vortices and streaks. Moehlis et al developed a 9-mode ODE 
model that was amenable to investigation in substantial detail. In particular they 
discussed the lifetimes of perturbations as the Reynolds number increased, locating 
the onset of complicated dynamics, and they also discussed the probabilistic distri- 
bution of lifetimes at fixed Reynolds number from randomised initial conditions of 
equal energy. 

It is clear that the assumptions made by these authors concerning spatial period- 
icity seems much easier to argue for in the streamwise and wall-normal directions 
than in the spanwise direction; indeed recent work by Schneider and co-workers 
(Schneider et al 2009; 2010a; 2010b) (see also Duguet, Schlatter and Henningson, 
2009) has concentrated on understanding the formation of structures which are 
spatially localized in z: very far from periodic in this direction. 

With this in mind we propose in this paper an extension of Waleffe's model 
which is a reduction of the Navier-Stokes equations to a collection of PDEs in z 
and t: we adopt the Fourier mode truncations that Waleffe used in order to remove 
the dependence on the streamwise x and wall-normal y coordinates since numerical 
work shows that, at least for some of these equilibrium and periodic orbit states, 
the flow structure in these coordinates can be well-approximated by a small number 
of Fourier modes. 

By retaining full dependence on the spanwise coordinate z we admit both the 
spatially-periodic solutions of Waleffe and the formation of localized states. In addi- 
tion, a wealth of spatio-temporal complexity is allowed. Our study is very similar in 
spirit to work by Manneville and co-authors (Manneville & Locher 2000; Manneville 
2004; Lagha & Manneville 2007) who preserve full resolution in two directions 
(streamwise and spanwise) and use Galerkin truncation only in the third ('wall- 
normal') direction. This enables these authors to consider the dynamics around 
turbulent spots that are localised in both the streamwise and spanwise directions, 
at the cost of more intensive numerical computations. Their work is therefore com- 
plementary to that presented here. We remark also that reduced models, of different 
kinds, have been used both in pipe flow (Willis & Kerswell 2009) and in under- 
standing the formation of turbulent-laminar bands in plane Couette flow (Barkley 
& Tuckerman 2007). 

The structure of the paper is as follows. In section [5] we discuss the derivation 
of our PDE extension of Waleffc's ODE model. In section |3] we present the results 
of our numerical investigations into the transition to turbulence described by the 
model. Wc conclude in section |4l 

2. Derivation of the PDE model 

In this section we define sinusoidal shear flow and summarise the Galerkin trunca- 
tion that we use to derive our simplified model for turbulent transition. 

Following Waleffe (1997) and Moehlis et al (2004), we use the usual Carte- 
sian coordinate conventions that x is the downstream direction ('streamwise'), y is 
the direction of the shear gradient, i.e. normal to the sidewalls and z is the span- 
wise direction. We write the Navier-Stokes equations for incompressible flow in the 
nondimensionalised form 

^ + u-Vu=-Vp + F(y) + — V^u, (2.1) 

ot Re 
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where we have scaled lengths by h/2 where h is the width of the channel, velocities 
by Uq the velocity of the laminar profile at a distance h/A from the upper boundary 
and pressure by U^p where p is the density of the fluid. The Reynolds number 
Re is therefore given by Re = Uoh/{2v), v being the kinematic viscosity. The 
evolution of (|2.1|) is subject to the usual incomprcssibility condition V • u = and 
the conditions for impermeable and stress-free upper and lower boundaries 



0, 



and 



dy 



duz 
dy 







at 



y = ±l. 



(2.2) 



The flow is assumed to be periodic in the x and z directions, with periodicities L^ 
and Lz respectively. We take the body force term F{y) to be 



F(2/) 



V2^ 
Re 



sinPye^, 



(2.3) 



where e^; is the unit vector in x direction and it is convenient to define the constant 
/3 = tt/2. The laminar profile is then the steady solution 



u = V2aml3yex, 



(2.4) 



of (|2.ip - (j2.3p . as shown in figure [TJ We note that although the 'sinusoidal shear 




Figure 1. Geometry of the domain, illustrating the basic sinusoidal shear profile 
u = ^/2 sin ISye^ which is sustained by the applied body force term. 



fiow' profile (|2.4p has an inflection point, the flow is linearly stable for all Re (Drazin 
& Rcid 1981). 

Wc now turn to our solution ansatz. We write the velocity field in the form 



um + V X ^T + V X V X $p, 



(2.5) 



where the subscripts denote mean, toroidal and poloidal components which are, 
respectively, expressed as sums of the first few Fourier modes in each case: 

Um = {Aisml3y + A2)e^, 
*T = ^3 cos I3y e^ 

+ (j44 sin ax — Ar, cos ax sin /3y — Aq cos ax + Ay sin ax sin /3j/) e^ , 
$P = As cos ax cos f3y By . 
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Table 1. Comparison of notation. 



This paper 


Ai 


A2 


^3 


A4 


^5 


^6 


Ar 


^8 


Waleffe (1997) 


M 


U 


V 


A 


C 


B 


D 


E 




const 


cos 72 


sin 72 


const 


const 


cos 72 


cos 72 


sin 72 



Wc define the streamwise and spanwisc wavcnumbcrs a ~ I-k jL^ and 7 = 2t: / L^ for 
notational convenience. The ampUtudes ^1, . . . , Ag are functions of z and t whose 
evohition can be obtained by substituting (|2.5p into (|2.1[) . Note that the form of 
the ansatz impUes that incompressibihty and the boundary conditions (j2.2p are 
automatically satisfied. 

The amplitudes Ai , . . . , Ag correspond exactly, in terms of their Fourier depen- 
dence in X and y, to the modes selected by Waleffe (1997). For reference, table [1] 
summarises the correspondence. 

We now briefly describe the role that each mode plays in the dynamics. Ai is the 
amplitude of the sinusoidal shear profile; the laminar state corresponds to Ai — \/2, 
A2 ^ ■ ■ ■ ^ Aq ~ Q. A2 describes variations in z of the streamwise velocity, i.e. 
the formation of streamwise streaks. A3 describes the formation of x-indepcndent 
streamwise vortices that redistribute the shear profile. Modes A4,. . . ,Aj describe, 
as in Waleffe (1997) the development of x-dependent distortions of the streaks, 
and in particular the linear instability of the x-independent streaks described by 
A2. These modes have no velocity component in the vertical (i.e. y) direction. Ag 
describes 'oblique rolls' and, in contrast to modes A4, . . . , Aj, has a non-zero vertical 
velocity component but no vertical vorticity component. 

To derive evolution equations for the modes Ai , . . . , Ag we use Fourier orthogo- 
nality in the x and y directions combined with projections onto individual compo- 
nents of either the velocity field u, the vorticity field a; = Vx u or (for Ag) the curl 
of the vorticity field. For later convenience we introduce the differential operators 
V^, Dp and P^o which correspond to the action of — V^ on different Fourier modes: 

We denote dAj/dz by A', for 7 = 1, . . . , 8. Considering first the e^ component of u 
in ()2.1[) we obtain the following PDEs in z and t for Ai and A2: 

dt + ^VJ^A, = -pA'2A, + ^iA'lA,~A4A'^) + ^{AeA'^~A'^Ar) 

+^(A'^A'g~a'A,A'g) + y^, (2.6) 

dt-^^d.?jA2 = -^{A,A,y + ^{A'lA,-Aa'^)-^{A'^A,-A,A'^) 

-^(A,Asy + ^iA',A'gy, (2.7) 

Now we turn to the vorticity uj. For A3, A4 and Aq we find that only one component 
of u) contains a contribution from each of these; taking the curl of (|2.5p the terms 
involving A3, A4 and Ag are 

u) = V p A^ cos I3y e^- + (V^A^ sin ax — V'^Aq cos ax) ey. 
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Therefore it is straightforward to consider the x and y components of the vorticity 
equation obtained by applying the operators e^; • Vx and e^ • Vx to (j2.ip in order 
to obtain evolution equations for Ay,, A4 and Ag: 



^Vl\vlAy = -o'piAiAj)' -a'PiA^Ae)' + ^i 



^ I jl^ A"\" , „ o2 A' M , ^ P 
2 



^{A,A'ir + a^'A'.A', + ^A,As 



■^A,Al (2.8) 



dt + ^VljVlA, = ^{A^A'i^A'(A,)~^A,A, + a{A,A'^-A'^A,) 

-a^A^Ae + ^{A^A'^)" - ^A^A', ~ a'pA'.Ar 

a^B'^ a/32 
^A,As + -^{A,A'^ - A'^As), (2.9) 

dt + ^Vl\vlAe = ^A,Ar-^iA,A'^-A'{Ar)+a^^A[As 

+^A^A',- ^{A^A'^Y' + o^'A^A, 

-aiA^A'l ~ A'^Ai) + ^{A.A',)" ~ a^/SAi^A^ 

-^M^^. (2.10) 

Evolution equations for A5, Aj and Ag are derived similarly, using (for A^ and 
A^) the projection operator ej, • Vx and (for A^) the projection operator e^ • Vx Vx 
applied to (|2.1I) . The resulting evolution equations are 

dt + ^Vl^VlAr, = a''AiAi~a{AiA'i~A'iAi)+a''A2Aj 

-aiA2A'^ - A'^Ar) - /3(A^A^)' 
+a'pA'^As - a^PiA.AeY + P{A,A'^)\ (2.11) 
1 
Re' 

+a{A2A'i - A'^A,) + PiA^A'iy 
-a^PiA^A^y, (2.12) 

1-Vl^^Vl^VlAs = -2a'p{A',A, + A'^A,)-2aP'A',A', 

+a{A'lAi)" - a{a^ + l3^)A':^Ai. (2.13) 



9t + -i^'Dlp)vlA7 = -a^AiAG+aiAiA'^^A'lAe)~a^A2A5 



Equations (j2.6p - (|2.13|) are a closed set of nonlinear PDEs in z and t which form a 
truncated model of the dynamics of sinusoidal shear flow. Crucially these equations 
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satisfy two consistency checks: the nonhnear terms in (|2.6p - (|2.13|) conserve energy 
and so reflect the conservative nature of the full u-Vu nonlinearity in (|2.ip . Secondly 
this model reduces to the model of Waleffe (1997) in the special case in which the 
amplitudes Ai, . . . ^ Ag are taken to be periodic in the z-direction, after applying 
appropriate projection onto orthogonal Fourier modes in z. We discuss each of these 
important issues in more detail in the following subsections. 

(a) Nonlinear terms and energy conservation 

In this subsection we show that the nonlinear terms in (|2.6p - (|2.13p do not 
contribute to the energy budget for the dynamics. This is a physically crucial prop- 
erty in constructing any reasonable reduced model for shear flows: the only energy 
source term must be the body force F(j/) that drives the laminar flow, and the only 
source of dissipation must be (linear) viscous diffusion. 

We define the (dimensionless) total kinetic energy of the flow to be 

E = - u-udxdydz, (2.14) 

where O = [0,7r/a] x [—1,1] x [0,^2] is (one half of ) the domain (in x,y,z coordi- 
nates) occupied by the fluid. We substitute our ansatz (|2.5p into (|2.14p and carry 
out the X and y integrals, noting that Fourier orthogonality enables us to remove 
all cross-terms except those involving A7 and As since they have the same Fourier 
dependencies; their contributions to u are 

(— A'j ain ax sin l3y \ / — a l3As sin. ax sin l3y 

, and ug = {a'^ As — A'^) cos ax cos Py 

a At cos ax sin Py J \ —fHA'sCosaxsinPy 

Hence the kinetic energy E is given by 

E = ^ T' Al + 2Al + {A'^f+p^Al + {A',f + a^Al 
^a Jo 

+ l{A',f + ^Al + iA',r + a'Al + \ {aPAs - Al.f 
+ 1 {a'As - Mif + \ [aA^ - {iN^f dz, 

= ^ / ' E{z,t)dz, (2.15) 

2a Jo 

which defines the 'local' energy quantity E{z,t). After integrating several terms by 
parts, and noting that the boundary contributions vanish (since we use a periodic 
boundary condition in z), and also after some manipulation of the cross-terms 
involving Ay and As, we obtain the expression 

^ = ^ / ' Al + 2Al + AsV}A3 + AiVlAi + ^A5VlA5 
+AeVlAe + -A^ViA^ + ]^AsVlpVlAs dz. 
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We can now compute the time evolution of the kinetic energy by difFerentiating 
with respect to time. After carrying out further integrations by parts we find 

JL = ^ j 2A1A1+4A2A2+2A3VIA3 + 2A4VIA4 
at 2a Jq 

+A5VIA5 + 2AgVIAg + AjVlA^ + A^VlplA^ dz. (2.16) 

We now substitute for the time derivatives using the PDEs (|2.6p - (|2.13p . The 
interest in pursuing this calculation is that at this stage we find that the conservative 
nature of the quadratic nonlinearities in (|2.6p - (j2.13|) becomes apparent since all 
the cubic terms cancel (after appropriate integrations by parts). We are then left 
with a single linear source term and a collection of quadratic dissipation terms: 

dE nP^V2 r^ A,dz~-^^ t 2 {^^ A\ + (Al)^) + 4{A',r 



dt aRe Jq 2aRe 

+2{'D}Asf + 2{VlA4f + {VlA^f + P^ {a'Al + {A',f) 
+2{VlA^f + {VlA^f + 13^ {o?A^, + (A'7)2) 
W^Vl^A^f ^{VlpJ^^f dz. 

The derivative of this equation for the evolution of the total kinetic energy, which 
describes the balance between the driving provided by the body force term F(y) and 
viscous dissipation, makes it clear that the nonlinear terms in (j2.6|) - (|2.13p conserve 
energy. This justifies the self-consistency of the selection of just the eight modes 
Ai, . . . , Ag in the reduced model: by including exactly these modes, no unphysical 
effects are introduced into the energy evolution. 

(6) Relation with the model of Waleffe (1997) 

Waleffe (1997) considered a far simpler modal truncation in which each mode 
contains only a single Fourier mode in z. This introduces an additional parameter: 
the wavenumber 7 in the z-direction, but the model comprises ODEs rather than 
PDEs for the eight amplitudes and is therefore much more amenable to analysis. 
The derivation of such a modal truncation is slightly simpler than that described 
above since one can use Fourier orthogonality directly in the z-direction as well 
as in X and y. In fact, the reduced model (|2.6p - (|2.13|) reduces exactly (after 
rescalings of the amplitude variables) to that discussed by Waleffe if one inserts the 
corresponding trigonometric dependencies and then projects out unwanted modes 
that arise from some of the nonlinear terms. This projection step is self-consistent in 
the sense that in both (|2.6p - (|2.13p and the ODEs derived by Waleffe the nonlinear 
terms conserve energy. Table [1] lists the Fourier dependence in the z-direction that 
Waleffe (1997) assumed for each mode. 

Waleffe observed that his ODE model had a number of deficiencies, for example 
the streaks (described by A2) were unstable to a;-dependent perturbations, with 
even and odd symmetry in y, only if the wavenumber 7 was sufficiently large: 
7^ > a"^ and 7^ > a'^ + /3^ , respectively, to be precise. 

More seriously, he discusses the inadequacy of his ODE model in capturing the 
interaction between the mean shear (M, or equivalently ^i) and the cc-dcpendcnt 
modes (A, C, B and D, or equivalently A4, . . . , Ai) that arise from consideration 
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of advection by the mean shear. We observe that in (|2.6p - (|2.13p these quadratic 
interactions take far more complex forms that in many cases vanish idcnticahy when 
the simple Fourier mode dependencies in z that are listed in table [T] are imposed. 
In particular the mean shear Ai is influenced by new combinations of cc-dependent 
modes which do not arise in the ODE model because the expressions AqAj — AqA^ 
and A'lAc^—AiA'l which are present in (|2.6p vanish identically in the ODE reduction. 
Terms with a structure identical to this (i.e. of the form AnA'^ — A'^Am) appear 
in several of the the other amplitude equations. Compared with Waleffe's ODE 
model, the other qualitatively new couplings introduced in (j2.6p - (|2.13p are the 
term PiAsA'^Y in ((2JT|) and the term -2a'^l3A[Ae in ((2l^ . 

3. Numerical results 

In this section we present the results of time-stepping the system of PDEs (|2.6p - 
(|2.13p over the range 50 < Re < 200. Our numerical method is the pseudospectral 
exponential time-stepping scheme referred to as 'ETD2' by Cox & Matthews (2002), 
using 128 Fourier modes in z, and suitably small timcstcps such that our results 
were insensitive to the timestep used. 

As in Moehlis et al (2004) our primary interest is in the emergence of a chaotic 
saddle in phase space. This is indicated by the lifetimes of chaotic transients as 
trajectories evolve towards the linearly stable laminar state, having started from 
initial conditions that are far from the laminar equilibrium. 

(a) Initial conditions and domain parameters 

The lifetimes of chaotic transients arc of course sensitive to the choice of initial 
condition, and in order to investigate the response of the flow to a spatially localised 
perturbation we took initial conditions corresponding to a Gaussian profile, scaled 
so that the four modes A3, A4, A^ and Aq gave equal contributions to the initial 
kinetic energy Eq. Later, in subsection [HIP, we compare these results with those 
obtained using a sinusoidal initial condition. 

Specifically our Gaussian initial condition takes the form Ai ~ A2 = Aj = Ag, = 
and 

^,=c,exp(-(i^||/^), (3.1) 

for J = 3, . . . , 6, where the coefficient a describes the width of the Gaussian, and 
the normalisation constants Cj are given by 



/ Epaa 



C5 = \ —1^7-T-^^T- (3-4) 



TTy/n (a^a^ 
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The differencics in the expressions for C3 , . . . , cg reflect the different contributions 
made by ^3, . . . , Ag to the kinetic energy (|2.15|) . 

We present results for a domain size L^ = I.TStt. Lz = 1.27r corresponding to 
the 'minimal flow unit' identified by previous authors (Hamilton, Kim & Waleffe, 
1995, and adopted by Moehlis et al 2004) as the smallest domain in which sustained 
spatiotemporally chaotic dynamics have been found numerically. We consider values 
for a in the range 0.2 < ct < 0.8, initially setting a = 0.2, so that the initial Gaussian 
disturbance is always spatially well-localised in the z direction. 

(&) Results at fixed Reynolds numbers 

Our numerical integrations are carried out until either the dynamics approaches 
very close to the laminar state, or until 1000 dimensionless time units of h/{2Uo) 
have elapsed: this is the maximum lifetime that transients are followed for in our 
computations. We follow transients for the range 50 < Re < 200, increasing Re in 
steps of unity, and increasing initial kinetic energy Eq in the range < E'o < 3.0 in 
steps of 0.01. 

At low Reynolds numbers. Re < 117, we find numerically that all initial condi- 
tions decay monotonically towards the stable laminar equilibrium, with a lifetime 
that increases slowly with Eq in the range < Eq < 1.7 approximately, and then de- 
creases slowly as Eq increases further. For Re> 117 the dynamics changes abruptly 
and initial energies around Eq = 1.7 show far longer transients, as indicated in fig- 
ure [D which shows the lifetimes of trajectories for four fixed values of Re, as Eq 
varies. The overall impression of figure [2] is of a well defined transition from laminar 
to spatio-temporally complex dynamics for a range of initial perturbations. It is 
clear that windows of rapid attraction to the laminar state remain, for example for 
Re ~ 143 and initial energies around Eq ~ 1.3. 

(c) Transition boundary and structure of the chaotic saddle 

The well defined 'transition boundary' at which there is an abrupt increase in 
lifetimes is shown clearly in figure [3] which presents a colour-coded surface plot 
of lifetimes as Re and Eq vary. In contrast to previous similar plots, for example 
figure 6 of Moehlis et al. (2004), and see also figure [TT] in the present paper, where 
spatio-temporally complicated transients appear at around Re = 120 at the ends 
of 'wispy' dendritic fingers that coalesce as Re increases, figure |3] indicates the 
sudden appearance of a 'fatter' chaotic saddle in phase space. The chaotic saddle 
appears to become less 'porous' as Re increases, and there appears to be only one 
substantial hole, at around Re = 140, Eq = 1.4. As in Moehlis et al. (2004) we 
might anticipate that figure [3] has fine scale structure as we zoom in. In contrast to 
the results of that paper we find however that finer-scale investigations do not reveal 
any coherent organisation to the lifetimes: instead of regions of concentric bands 
of different colours, for example, we find merely fine-scale apparent randomness 
and intermittency. This is illustrated in figure IH Figure EJ^a) shows the region 
130 < Re < 140, 1.6 < Eq < 1.7 with a ten-fold increase in the resolution along 
each axis. Figure|Ub) shows a further order of magnitude increase in the resolution 
both in Re and in Eq; this part of the figure corresponds to the region within 
the black box in the centre of figure Hl[a) . This very detailed enlargement appears 
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Figure 2. (Online version in colour.) Lifetimes of transients started from the Gaussian 
initial condition, for initial kinetic energies < -Eo < 3.0 for Re — 108 (black, +), 
Re =117 (red, o), Re = 128 (blue, 0) and Re = 143 (magenta, D). Parameters are 
a = 0.2, L^ = 1.757r, L^ = 1.27r. 



to exhibit some slight correlation, for example a propensity to favour lifetimes of 
around 600 time units at Re « 135.4 but otherwise the unpredictable intermittent 
nature of the dynamics continues as we consider finer and finer divisions in Re 
and Eq ■ This qualitative observation should be contrasted with the lifetimes figures 
plotted by Moehlis et al. (2004) which exhibit much more structure. We return to 
this point in section [21|f| . 

The lower boundary of the chaotic saddle is of particular physical interest, since 
it describes the smallest amplitude perturbation required to produce an extended 
chaotic transient. For the particular form of perturbation used here, in the case 
a = 0.2, we find numerically that the lower boundary in figure[3]is remarkably well 
fitted by a power law curve as illustrated in the log- log plot in figure [5] Additional 
computations show that a power law continues to provide a very good fit as i?e 
is increased, up to at least Re — 10"^, as shown in figure [BJ^a). The black curve in 
figure [BJ a) is the best-fit curve 



Sn = 



122 



4.6 



(3.5) 



Blue squares in figure IHl indicate that the flow returned to laminar before the com- 
putational time limit T,„ax was reached, while red + signs indicate that the flow 
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Figure 3. Lifetimes of transients started from the Gaussian initial condition, over the 
range Q < Eo < 3.0 and 50 < i?e < 200. Parameters are a = 0.2, L^, = I.TStt, L^ = \.2-n. 



did not rclaminarise before t = Trr 



We set T,r 



5000 for 200 < Re < 750 and 



Trr 



10000 for 800 < Re < 950. 



In figure inijb) we show the best-fit power-law for the boundary between laminar 
and spatio-temporally complicated flow for different widths of the Gaussian initial 
condition p.ip . All are excellent fits to the data over the range 200 < Re < 1000, 
and the results are identical for computations using 128 modes or 32 modes in z. 
Since the slopes of the power-laws become more negative as the coefficient increases, 
the overall envelope of perturbation energies as Re increases, formed by taking the 
minimum value over all four curves, has contributions from each curve at sufficiently 
large Re. This envelope is an estimate of the true lower boundary for which we would 
ideally compute the minimum over all possible forms of initial perturbation. 

For perturbations with a Gaussian profile p.ip , and for the range of a considered 
here, it appears unlikely that the differences in exponents would be detectable 
experimentally over the range of relevant and accessible Re. Table [5] gives the 
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(a) (b) 

Figure 4. (Online version in colour.) Enlargements of figure [3] showing a lack of coherent 
organisation to the lifetimes at higher resolution in Re and Eq. (a) Transient lifetimes 
over the range 1.6 < Eo < 1.7 and 130 < Re < 140. (b) Transient lifetimes over the range 
1.65 < Eo < 1.66 and 135 < Re < 136 as indicated by the black square in (a). For both 
figures, the other parameter values are a = 0.2, Lx ~ 1.757r, L^ = 1.27r. 

details of the best-fit parameters of the power-laws shown in figure IH^b) , using a 
functional form Eq = [c/ReY . 



(d) Lifetime distributions 

Moehlis et al also computed the distribution of lifetimes of turbulent transients 
at fixed initial energies and Reynolds numbers. They found that the survival proba- 
bility P{T) that the solution had not decayed back to the laminar state after a time 
T was distributed exponentially, with a mean lifetime that increased rapidly with 
Re above the critical value at which the chaotic saddle appeared. Numerically, we 
construct a lifetime distribution by distributing the initial energy randomly across 
a subset of the modes. For our spatially-extended model there are two possible 
approaches to randomising the distribution of energy. 

In the first approach, the total initial energy is distributed uniformly on a spher- 
ical energy shell. This can be achieved straightforwardly by modifying the expres- 
sions (j3.2[) - p.4p for the normalisation coefficients by replacing Eq with 4£'o^| in 
the expression for coefficient Cj, where j = 3, . . . , 6. The coordinates ^j are those of 
points distributed at random over the unit sphere in M^ given by ^j=3 Cj = 1- We 
refer to this as randomising over the amplitudes of the perturbations. 

In the second approach, the central position of each Gaussian perturbation Aj 
is shifted randomly in z away from the centre of the domain Lz/2, i.e. the values of 
the coefficients C3, . . . , cg are left unchanged but the centre of the Gaussian in p.ip 
is modified for each amplitude Aj. Since we employ periodic boundary conditions, 
the total energy in the perturbation is preserved. We refer to this procedure as 
randomising over the locations of the perturbations. Figure [7] shows the results of 
the first approach at two Reynolds numbers quite close to the transition boundary. 
At low lifetimes (T < 200) the survival probability remains unity since any transient 
takes at least this finite amount of time to decay to within the prescribed threshold 
of the laminar state. As T increases there is an initial sharp drop indicating that a 
substantial proportion of trajectories do decay rapidly, with lifetimes in the range 
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Figure 5. Log-log plot of transient lifetimes over the range 0.1 < E < 3.0 and 
100 < Re < 200 to show power law form of the lower boundary of the attractor. The 
solid (red) line indicates the power-law Eq = (^) . Parameter values are a — 0.2, 
La: = 1.757r, Lz = 1.27r. 



200 < T < 400. At larger T the distribution becomes close to a straight line on 
the linear-log plot, which is consistent with an exponential distribution P{T) ^ 
ai cxp(— a2T) for T > 400 (approximately), as indicated by the dashed line. For 
Re — 130 the best-fit coefficient values are ai = 3.93 x 10~^ and 02 = 8.65 x 10~^. 
For Re = 140 the best-fit coefficient values are ai = 4.57 x 10~^, 02 = 5.11 x 
10"''. Similarly, figure |S] shows the lifetime distribution computed from the second 
approach, randomising over the locations of the perturbations. The distribution 
shows very similar characteristics to that in figure UI&) and is again well-described 
by an exponential distribution with best-fit parameters ai = 4.24 x 10~^ and 02 = 
8.32 X 10~^. We observe that the values for the coefficient 02 are very similar 
between the two randomisation methods. In all cases, we expect that as Re increases 
the distribution exhibits systematic deviation from exponential and flattens out. 
We anticipate that, as in the ODE case, this is due to the appearance of stable 
attractors near the chaotic saddle, as Re increases. These attracting sets then absorb 
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(a) (b) 

Figure 6. (Online version in colour.) Estimates of the lower boundary of the chaotic saddle 
over the range 200 < Re < 1000. (a) a = 0.2; squares (blue) indicate rapid return to the 
laminar state. '+' symbols denote a long-lived chaotic transient. The black line indicates 
the power-law scaling of the lower boundary, (b) Best-fit power laws for a range of values of 
a, showing a broad insensitivity to the width parameter. Parameter values are L^ = 1.757r, 
L, = 1.27r. 



Table 2. Dependence of power-law scalings on the width a of initial condition, 
a Constant, c Exponent, p 
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Figure 7. (Online version in colour.) Distribution of lifetimes P{T) computed over an 
ensemble of 2000 initial conditions produced by randomising over amplitudes of perturba- 
tions at Eo = 1.0. (a) Re = 130 (b) Re = 140. Parameter values are a = 0.2, L^ = 1.757r, 
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Figure 8. (Online version in colour.) Distribution of lifetimes P{T) computed over an 
ensemble of 2000 initial conditions produced by randomising over the locations of pertur- 
bations at Eq = 1.0 and Re = 130. Parameter values are a = 0.2, Lx = 1.757r, Lz = 1.27r. 



trajectories with positive probability, leading to a proportion of trajectories which 
never return to the laminar state (Moehlis et al 2005). 



(e) Spanwise resolution 

The very high resolution used in the z direction {N — 128 Fourier modes) is 
greater than required, for such a small domain, in order to capture the dynamics 
of the 'active' modes of the system. In order to probe the range of wavenumbcrs 
that make substantial contributions to the dynamics we investigated the effect of 
varying the truncation level of the numerical scheme in z. We define the 'm-mode 
truncation' of the PDEs (|2.6p - (|2.13p by keeping the Fourier modes ~ e^™'''^ for 
< n < 771 — 2 for Al , ^4 and A5 , and keeping the modes for < n < tti — 1 for 
A2, A3, Aq, A^ and A%. The m-mode truncation is thus a collection of (at most) 
16m— 14 real ODEs. This is an upper bound on the effective dimension of the ODE 
dynamics since not every real and imaginary part is coupled for every m. 

By construction, this further reduction resembles the Waleffe model in the par- 
ticular case 771 = 2, but by varying m we are able to probe the influence of the 
higher-wavenumber modes on the location of the lower boundary for the onset of 
temporally complex dynamics. For each value of Eq and tti, a single initial condition 
was used. This initial condition was the projection of the Gaussian profiles described 
in section [3l^a) onto the available Fourier modes in z. Results are shown in figure 13 
Figureinia), for the small domain Lz = 1.2tt, shows that computations with trunca- 
tion levels m > 5 produce an identical indication of the boundary between laminar 
and spatio-temporally complex behaviour. The case m = 4 is qualitatively but not 
quantitatively correct. For the larger domain L^ = 67r, figure [9jb) indicates that 
the numerical results are essentially unchanged for m > 17. The number of Fourier 
modes required to maintain accuracy in the larger domain is therefore broadly in 
line with, although slightly lower than, what one might naively expect. 

Finally we note that as m decreases further, the boundary appears consistently 
to move to lower Eq. This behaviour is purely a function of the organisation of 
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Figure 9. (Online version in colour.) Lifetimes of transients T" as a function of initial energy 
-Bo for Ke = 200 in domains of widths (a) L^ — 1.2-n and (b) Lz = Gtt, showing the effect 
of varying the truncation level m. Vertical axis is logarithmic for clarity, and shows the 
additional lifetime after subtracting a constant: T — 365 in (a), T — 414 in (b). In (a) the 
data for m = 9 and tti = 5 lie exactly on the m = 16 values: they are artificially offset 
vertically for clarity. In (b) the data for m = 19, m = 17 and in = 15 lie exactly on 
the m = 32 values and are similarly offset for clarity. Computations were terminated at 
Tmax = 5000. Parameter values are a = 0.8, L^ = 1.757r. 



invariant sets in phase space; it is not clear that there is a physical reason why 
these should appear to move in one direction or the other. 

(/) Spatio-temporal dynamics and effect of initial conditions 

In this final subsection we comment on the spatio-temporal dynamics of the 
PDEs for initial energies near the transition boundary, and we compare the results 
of section [3l^c) for the lifetimes of transients with those in this section obtained us- 
ing sinsoidal initial conditions. Figure [TU] shows the spatial and temporal evolution 
of the PDEs for Re = 200 and £'o — 0.11, just above the transition boundary. The 
quantity E{z^ t)-- A\ + {Ai — -\/2)^ is plotted in the figure, so that the laminar state 
Ai ~ \/2 is at level zero, where E{z,t) is defined in (|2.15p . The initial monotonic 
decay towards the laminar state is interrupted by the growth of oscillatory distur- 
bances. These disturbances generate a sharp spike in the kinetic energy, localised 
both in space and time, before the solution settles into a spatio-temporally chaotic 
state. For comparison, at Eq = 0.1 we observe only monotonic decay to the laminar 
state. 

Analysis of the energy distribution across Fourier modes shows that the small- 
scale oscillations just prior to the localised spike involve many higher-wavenumber 
Fourier modes. These mode amplitudes grow very rapidly as we approach the spike 
state. Since the formation of the spike is, in almost every case, the 'edge state' that 
is the precursor to spatio-temporally complicated dynamics, one interpretation of 
the results in the previous sub-section is that the spatial resolution required to 
correctly determine the boundary of the basin of attraction of the laminar state is 
indicated by the resolution required to describe accurately these spatially-localised 
spikes. 
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Figure 10. (Online version in colour.) Space-time plot of the local energy quantity 
E{z, t)—A\ + {Ai — a/2)^, indicated by both surface height and colour, for the PDEs (|2.6|l - 
(|2.13|l for Re — 200 in the minimal flow unit domain. E{z, t) is defined in equation (|2.15|) . 
The additional square terms shift the laminar state to the zero level in the plot and serve 
to highlight the localised 'edge' state at t ~ 450. Parameter values are Eo — 0.11, a — 0.2, 
L^ = 1.757r, L, = 1.27r. 



More generally, the dynamics near the boundary between monotonic relami- 
narisation and spatio-temporal complexity appear to depend on the evolution of 
high-wavenumber modes, seeded by the use of a Gaussian initial condition which 
injects energy into every available mode. This is in contrast with the sinusoidal 
initial conditions used in previous reduced models, where, by construction, such a 
sinusoidal initial condition was the only possible choice. 

In figure [TT] we show the analogous plot to figure [3] for the lifetimes of transients, 
but in this case using a sinusoidal initial condition instead of a Gaussian profile. 
For the sinusoidal initial condition we set A2 = C2COS7Z; A3 = C3sin7z; A4 = C4; 
A^ = c^; Ai ~ Aq = Aj = Ag, = where the constants C2 , . . . , C5 are given by 



C2 



aEo 
27rL, 



aEn 



C3 



(/32 -f J^)7tL, 



C4 



E, 







2TTaL, 



C5 



En 



airLz 



so that the inital kinetic energy £^01 defined in (|2.15|) . is distributed equally between 
the four non-zero amplitudes. 

Comparing figures [3] and [TT] there are important qualitative and quantitative 
differences. Firstly, the boundary between relaminarisation and spatio-temporal 
complexity is much more obvious in figure [3] where we employ the Gaussian ini- 
tial condition. The boundary in figure [TT] shows much more 'structure'; it is much 
less obvious that a boundary in the sense, for example, of a countable collection 
of (piecewise) continuous curves in the {Re,Eo) plane, can even be defined. Many 



Article submitted to Royal Society 



Transition in a ID shear flow model 



19 



Ui 1.5 




1000 



900 



800 



700 



600 



500 



400 



300 



200 



100 



Figure 11. Lifetimes of transients started from the sinusoidal initial condition, over the 
range < -Eo < 3.0 and 100 < Re < 200. Solid (red) line indicates the boundary 
-Bo — (75.7/-Re)^''* corresponding to the lowest of the boundaries in figure [6l^b) for Gaus- 
sian perturbations. Parameter values are L^ = 1.757r, L^ — 1.2ii. Computations used 
A^ = 32 Fourier modes in z and were terminated at Tmax = 1000. 



deep valleys of short lifetimes persist up to Re ~ 200 and beyond. Secondly, the 
solid (red) line in figure [TT] shows the lowest boundary computed for Gaussian per- 
turbations, corresponding to cr = 0.5, see table [2] It appears that that the laminar 
state is substantially more sensitive to Gaussian perturbations than those of the 
same energy but in the form of low-wavcnumber sinusoids. Tentatively, based only 
on the data in figure 1111 wc suggest that the lower boundary of the appearance of 
spatiotcmporally complex dynamics arising from the sinusoidal perturbation scales 
as Eq ^ Re^^ , i.e. perturbation amplitude scaling as Re~^ . This is a typical expo- 
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nent produced by very many reduced models of very low order, as summarised and 
discussed by Baggett & Trefethen (1997). 

Within the region of spatio-temporally complex dynamics, the lifetimes in fig- 
ure [TT] and in additional enlargements (not shown here) show a degree of correlation 
between lifetimes at neighbouring points in the (i?e,i?o) plane which is not nearly 
so clearly shown to exist in figure [3] or in the enlargements shown in figure SI 

In summary, in these respects figure [TT] is reminiscent of lifetime plots for the 
very low-dimensional models of Moehlis et al (2004), see their figures 6 and 8, 
and figures 5 and 6 in the paper by Eckhardt & Mersmann (1999). We conclude 
that allowing for the accurate representation of spatially-localised initial conditions 
by extending the spanwise resolution of the model generates results that differ 
substantially, both qualitatively and quantitatively, from those of previous reduced 
models. 

4. Discussion and conclusions 

In this paper we have presented an extended version of the Galerkin-truncated 
model due to Waleffe (1997) for the transition to turbulence (or, at least, spatio- 
temporally complex dynamics) in sinusoidal shear flow. This model is appealing 
since it provides an intermediate step between previous analytical work and DNS 
of the full Navier-Stokes equations. Preserving full resolution in the spanwise {z) di- 
rection and removing the assumption of periodicity allows both the use of spatially- 
localised initial conditions, and the (transient) formation of localised structures in 
the flow which (although unstable) are known to exist and play a role in mediat- 
ing the onset of turbulence. The use of a small number of Fourier modes in the 
wall-normal (y) and streamwise (a;) directions provides the simplification of the un- 
derlying Navier-Stokes equations, which in turn allows us to perform very detailed 
investigations of the dynamics of this reduced model. 

We compare our results with those of previous authors, in order to see which 
properties are common to these different approaches, and which are not. For exam- 
ple, we find, in agreement with the results of Moehlis et al 2004, that the lifetimes of 
turbulent transients are well-described by an exponential distribution. However, our 
results show that the transition boundary, while exhibiting some of the 'structured' 
shape observed by many authors (including, in the case of pipe flow, Schneider, 
Eckhardt & Yorkc 2007), appears at lower perturbation energies, and much more 
abruptly, than for the ODE models investigated by Eckhardt & Mersmann (1999) 
and Moehlis et al (2004) . The PDF model that we present here is able to represent 
both spatially-localised and spatially-extended initial conditions and therefore we 
are able to make direct comparisons of this kind. 

Our key finding is that spatially-localised initial conditions are able to provoke 
complicated behaviour at substantially lower energies than the sinusoidal, spatially- 
extended perturbations used in previous studies. Moreover, the perturbation energy 
at the lower boundary of the chaotic saddle appears to scale as Re^ with the expo- 
nent p ~ —4.3, rather than Re~^ as in Eckhardt & Mersmann's 19-mode truncated 
ODE model (note that their figure 5 showing an Re~^ power law plots Re against 
mode amplitude which is proportional to E'q ). In the present work, the exponent 
in this power-law scaling was found to depend only weakly on the width of the 
Gaussian perturbation used. 
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In addition, our results are robust to the numerical resolution used in the span- 
wise direction, and, for a relatively small domain of width Lz = 1.27r, point to the 
necessity of keeping around 5 Fourier modes in z in order accurately to capture 
the dynamics of the fully-resolved PDE model. One possible explanation of these 
results is that admitting higher-wavenumber modes generates many more invariant 
sets within the boundary of the basin of attraction of the laminar state. Then, even 
small amounts of initial energy in these modes forces the system to spend much 
longer in the vicinity of these sets before being able to relaminarise. In this sense, 
'holes' in the basin boundary are filled in. The existence of these new invariant 
sets, and the lack of 'holes', leads to a robustness in the lengths of transients, and 
therefore to a more clearly defined boundary between monotonic rclaminarisation 
and longer-lived transients. 

It would clearly be of interest in future work to look at the relation between 
localised states which have been observed and studied in some detail in DNS for 
shear flow problems (Schneider et al 2010a, 2010b) and the dynamics of the reduced 
model presented here. We anticipate that the reduced model contains such states, 
and the homoclinic snaking bifurcation diagrams that typically organise them in 
driven dissipative systems such as shear flows, just as model ODE truncations, for 
example that discussed in Moehlis et al 2005, contain equilibria and time-periodic 
solutions very similar to those located in DNS (Nagata 1990; Gibson et al 2009). It 
should be possible systematically to further reduce the model equations presented 
here in order to make direct connections between theoretical work on localised 
states (Burke & Knobloch 2006; Chapman & Kozyrcff 2009; Dawes 2010) and the 
DNS results referred to above. In turn, the identification and analysis of additional 
unstable invariant sets within the boundary of the basin of attraction of the laminar 
state (as discussed by Lebovitz 2009), and their parameter dependence, will greatly 
help our understanding of the process of relaminarisation. 

JHPD would like to thank Rich Kerswell and Tobias Schneider for useful conversations, 
and Matthew Chantry for a minor correction. Both authors are grateful to the anonymous 
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